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SUMMARY 

An automatic version of the multigrid method for the solution of linear systems arising from the 
discretization of elliptic PDE’s is presented. This version is based on the structure of the algebraic 
system solely, and does not use the original partial differential operator. Numerical experiments show 
that for the Poisson equation the rate of convergence of our method is equal to that of classical 
multigrid methods. Moreover, the method is robust in the sense that its high rate of convergence is 
conserved for other classes of problems: non-symmetric, hyperbolic (even with closed characteristics) 
and problems on non-uniform grids. No double discretization or special treatment of sub-domains 
(e.g. boundaries) is needed. When supplemented with a vector extrapolation method, high rates 
of convergence are achieved also for anisotropic and discontinuous problems and also for indefinite 
Helmholtz equations. A new double discretization strategy is proposed for finite and spectral element 
schemes and is found better than known strategies. 


INTRODUCTION 


The multigrid method is a powerful tool for the solution of linear systems which arise from elliptic 
PDE’s [1] [2]. This is an iterative method, in which the equation is first relaxed on the original fine 
grid in order to smooth the error; then the residual equation is sent to a coarser grid to be solved 
there and to supply a correction term. Recursion is used to solve the coarser grid problem in the 
same way the original equation is handled. In order to apply this procedure, the differential operator 
has to be discretized on all grids, and restriction and prolongation operators have to be defined 
in order to pass from coarse to fine grids and vice versa. The multigrid method works well for the 
Poisson equation in the square, but difficulties arise with non-symmetric problems, indefinite problems 
and problems with discontinuous coefficients or non-uniform grids. In those cases, it is not easy to 
discretize the differential operator on coarse grids and to generate the restriction and prolongation 
operators appropriately. Some suggestions about how to handle discontinuous coefficients are given 
in [4] and [5], while the singularly perturbed case is discussed in [6]. Slightly indefinite problems are 
discussed in [7]. These approaches involve special treatment of problems according to the original 
PDE, and the need for a uniform approach is not yet fulfilled. 
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In principle, the multigrid procedure is problem-dependent, and cannot serve as a black box 
that solves every problem. Special attention has to be given to the neighborhood of the boundary 
and to lines of discontinuity. In [8] [9] [10] an algebraic multigrid method for symmetric problems 
is developed. Though this method is automatic in the sense that it depends on the linear system 
of equations solely, it suffers from the disadvantage of the coefficient matrices for coarse grids being 
of 9-diagonal type, even when the original matrix is of 5-diagonal type [11]. An algebraic version 
of multigrid which overcomes this difficulty is presented in [12], and generalized to nonsymmetnc 
problems in [13]. This version, however, does not improve the classical multigrid in cases of indefinite 
or hyperbolic problems and of non-uniform grids. 

The algorithm which is presented in this work, and which we denote Multi Block Factorization 
(MBF) (the reason for this terminology will become clear in the next section), gives a uniform 
approach that enables one to handle tne above difficulties. It relies on the algebraic system, of 
equations solely, and not on the original PDE. The operators for coarse grids, as well as restriction 
and prolongation operators, are automatically defined when the coefficient matrix is given. It seems 
to be more robust than the classical multigrid method, as it solves non-symmetric problems (even 
hyperbolic or with closed characteristics) as quickly as classical multigrid solves the Poisson equation. 
Moreover, it is applicable to non-uniform grids as well, and does not require any special treatment 
of sub-domains. For anisotropic, discontinuous or indefinite problems M BF by itself is not always 
sufficient. However, it can cope with such problems successfully when it is applied in conjunction 
with vector extrapolation methods. In our numerical examples in the present work we have employed 
the Reduced Rank Extrapolation (RRE) of [17] and [18]. This and other related methods have been 
surveyed in [19] and their analysis provided in [20], [21] and [22]. The numerical implementation that 
we have used is the one given in [23] . 

The MBF algorithm is described in Section 2. In Section 3 numerical results are presented. In 
Section 4 the algorithm and the numerical results are discussed. 


2 DESCRIPTION OF ALGORITHMS 
2.1 Definition of the TBF Method 

Let A be an N x N matrix. Let x and b be iV-dimensional vectors. Consider the problem 

Ax = b (1) 

An iteration of the Two-Block Factorization (TBF) method is defined by 

TBF(x in , A, b, x out ) : 

20 = 

x <+ i = Xi - Zi(Ax{ - b ) 0 < i < i\ 

Qe = R(Axi 1 - b) (2) 

x il+ i = Xj, - Pe 

Xj+i = Xj — Zi(Axi - b) ii < i < ii + k 

%QUt ~ + 

where the Z { are some preconditioning operators, k and k are nonnegative integers denoting the 
number of presmoothings and post-smoothings respectively and R , P and Q are operators to be 
defined later. Define __ 

® out = *^out 


■ w !•*.- 
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Then 


»l + «2 *» — 1 

«.«*=[ n (I - ZiA)](I - PQ-'RAXUil - ZiA)\e in 

»=ti+l i=0 

Consequently, 

Q = RAP =4* e ou t = 0 =>■ x ou t = x. (3) 

In the sequel, practical choices for Q will be considered. An iterative application of TBF is given 

by 

xo = 0, i = 0 

while ||rest<fua/|| > e 

TBF(xi, A , 6, x<+i) 
i *— i + 1 

end while 

2.2 Definition of the MBF Method 

The Multi-Block Factorization (MBF) method is a modification of the TBF method, in which the 
system (2) is not solved directly, but is divided into several independent subsystems, which are solved 
directly or recursively by MBF itself. For simplicity, we first write the algorithm for tridiagonal 
systems. The operators P, R , Q and D will be defined later. 

MBF(x in , A, 6, x out ) : 
if A is diagonal 

Xout = A~ l b 
otherwise: 

D = diag(di,. . . 

MBF (0, D~ x Q , D' 1 R(Axi n — b), e) 

Xout — X{ n Pe. 

We turn now to the more general definition of MBF. First we note that if there exists a subset 
of coordinates of x which are independent of the others, then there exists a projection II onto the 
sub-space spanned by those coordinates such that (nAtl)IIx = 116. In the following definition of 
MBF such sub-systems are solved directly, provided this can be done easily. The co-subsystem is 
solved recursively. 

MBF(x in , A, 6, Xout) : 

1. If A is diagonal, set x ou t = A -1 6 and stop. 

2. If A includes an independent tridiagonal subsystem (IIAnjllx = 116, solve it directly: IIx out = 
(IIAn) -1 n&. If not, set II = 0. 

3. 

yo = (I - n)x in 
6 = (I-U)b 

J/i +1 = Vi -(I- n )Zi(Ayi - 6) 0 < i < i\ 

D = diag(d\, . . . ,dff) 

MBF (0, D~ l Q, D~ x R(Ay il - 6),e) (4) 
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y.i+1 = y.i - (i - n)Pe 
y,+i = yi-(I - n)Zj(.Ayi — fc) *1 < * ^ *i *b *2 

(/ n)x 0U f = Vh+ii+i- 

Trivially, one can replace action (4) by variant a : 

MBF (0, QD~\R(Ay h - b), e) 
e = D~ l t 

or variant 6 : 

MBF (0,r>- 1/2 <5r>- 1/2 ,^ _l/2 i?(Ayn - 6),e) 
e = 2T 1/2 e. 

An iterative application of MBF is given by 

xo = 0, i = 0 

while ||residua/|| > e 

MBF{xi,A,b , x,+i) 
i <— i + 1 

end while 

2.3 The Tridiagonal Case 


Let I denote an identity operator. Suppose N = 2 n for some positive integer n, and let A be a 
tridiagonal M-matrix satisfying diag(A) = I. Let M(N) be the permutation matrix which reorders 
the variables of JV-dimensional vectors such that odd numbered variables appear in a first block and 
even numbered variables appear in a second block. Define 

M 0 = M{N), A 0 = A. 


I -Bo 
0 I 


Then for some bidiagonal matrices Bq and Co we have 

Ai = ^ B f ) M„ = R- a ',iQa,<P& 

where . . 

= ( - Co / ) M °’ QaA = ( 0 I- Co Bo ) ’ Pa ' 1 = M ° T 
Note that Qa , i is the Schur complement for A. 

For i = 1,2,..., let i* denote an identity operator of order N — 2 n_ '. Let Mi be the N x N 
permutation matrix that reorders the coordinates Xj, i = N — 2 n_l + 1, . . . } N, of an A - - dimensional 
vector in the above manner, that is, order odd coordinates in a first block, then even coordinates in 
a second block. In fact, 

Mi = ( 0 M(2 n -') ) ' 
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For 1 < i < n, define 

Ii 0 0 \ / J, 0 0 \ 

Pam i = M? 0 I -Bi I , R Ati+l = 0 I 0 

0 0 I ) V 0 -Ci I ) 

( Ii 0 0 

Qam i = P-amiA\P ami = I 0 I 0 

V o 0 I- CiBi 

Dami = diag(Q A M 

( Ii + 1 0 0 \ 

Aj+i = Dx)i + iQA,t+i = -Wi+i 0 I B i+ 1 M i+ i, (i < n - 1). 

V 0 c i+ 1 / / 

The last equality sign implicitly defines B{+ \ and Ci+\. For variants a and b, the above definition is 
modified to read 

■Ai+i = Qami d ami 

and 

■At+i = dJ + 1 q AiM Da,U 

respectively. 

Lemma 1 For all variants , the matrices Qi and are tridiagonal M -matrices. 

Proof: The lemma follows from the definition by induction on i. □ 

Lemma 2 The TBF method, when applied with 

Q = Qa,u P — Pa, ii P = Ra.u *i = 0, *2 = 0 

is a direct method. 

Proof: Since 

Q = Qa, i = Ra,\MPa,\ — RAP 
the lemma follows from equation (3). □ 

The even numbered variables may be viewed as coarse-grid points. Then Q is a coarse grid 
operator, R is a fine- to- coarse restriction and P is a coarse-to-fine prolongation. 

Lemma 3 The MBF method applied with the operators 

A = Aj_i, Q = Q a ,„ D = D A< i, P = Pa,« R = Ra,{ 

on the i th call to the MBF procedure, is a direct method. 

Proof: Note that on the (n+l) 4< call to the MBF procedure, the coefficient matrix A n is diagonal, so 
the MBF procedure is a direct solve. By induction on i = n , . . . , 1, all calls to MBF are equivalent 
to calls to TBF , hence are direct solves. □ 

In fact, in the tridiagonal case the MBF method is equivalent to the cyclic reduction method. 
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Note that if the matrices P and R are defined to be the rectangular matrices 

Pam i = Mf ( = ( -C, l)M„ 

Qa,i+ 1 = RA,i+\AiPA,i+l — I ~ CiBi, 

Da ,\+ i — diag(QA,i+\), 

A,+i = -DlJi+iQx.i+i = ^ ^ B 'j 1 M t+ i, 

then the algorithm will still be applicable. As a matter of fact, the only difference between this method 
and the former is that in the present method we are not taking advantage of the known residuals on 
the odd numbered variables when making the coarse grid correction. Hence, if those residuals are 
zero, which may happen as a result of red-black presmoothing, the present method serves as a direct 
method, just as the former. 


2.4 The Separable 2-Dimensional Case 


Let S = (sjj) and T = (tij) be matrices of order M and N respectively. Define 


S oT = 


si t iT 


s \,mT 


Si,jT 


sm,iT 


sm,mT 


Actually, o denotes the tensor product. Suppose A is of the form 

A — T o E s ,o A Ej t o 0 S 

where T and S are tridiagonal scaled M-matrices and Es , o and Et,o are diagonal matrices. For 
example, if 

T = S = tridiag(- 1/2, 1, -1/2), E T , o = E s , o = / 

then A represents a central discretization of the Poisson equation on a square with Dirichlet boundary 
conditions. 

Suppose T and S have the same order N, which is a power of 2. As in the previous section, we define 
the matrices T*, Si, Rt,» Rs,h Pr,i, Ps,% , Qtj, Qs,i, Dx,i and Ds,»- For any matrix B = (oi,j)i<i,j<N 

^ N 

row3um(B ) = diag(%2 bi,j)i<i<N- 
i=i 
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For 0 < i < n, define 


Et = rowsum(RT,i+i)ET,i’ rowsum(PT,i+\) 

E s = rowsum(Rs,i+i)Es,i- rowsum(P Sl i+i) 

Pa,m = -Pr.i+i 0 -Ps.i+i 

Ra i»+i = -Rr,«+i o Rs,i+ 1 

D a i»+i = -Dr,»+i o Ds,i+ i 

J5 r ,t+i = ^rl+i-^r 

£?5,i+l = -C>S, 1+1^5 

Qx.t+i = -Rr.i+i^iPr.i+i o £5 + o Ps.t+iS'iPs.i+i 
A«+i = D^+xQaa+i 

— TJ+1 0 ■E'S.i+l + Ej,i+l 0 ‘S'i+1- 

In the definition of TBF, we take 

Q — Qa,u P — Pa,u R — P/1,1- 

Note that if we would eliminate the word rowsum in the above definition, the TBF method would 
be direct, due to equation (3). Nevertheless, for smooth vectors, multiplication by a positive matrix 
B is well-approximated by multiplication by rowsum (B). Since TJ is an M-matrix, Rr,i+i and Pj il+ i 
are positive. Consequently, if we use presmoothing and post-smoothing, i.e. ii > 0 and it > 0, 
then the error is smooth, so equation (2) would give a good corrector for the current approximation. 
Moreover, the use of the above row-sum approximation makes the system (2) much easier to solve 
than the original system, since it includes 4 independent subsystems: 

1. A diagonal system connecting variables which are odd in both directions, i.e., variables that 
correspond to odd rows of both S and T (fine grid system). 

2. A tridiagonal system connecting variables which are odd in the first direction and even in the 
second, i.e., variables that correspond to odd rows of T and even rows of S (half coarse system). 

3. A tridiagonal system connecting variables which are even in the first direction and odd in the 
second, i.e., variables that correspond to even rows of T and odd rows of S (half coarse system). 

4. A penta-diagonal system connecting variables which are even in both directions, i.e., variables 
that correspond to even rows of both S and T (coarse grid system). 

Only the solution of the last subsystem is expensive. The MBF method solves this subsystem 
recursively by the same procedure. On the i th call to MBF (0 < i < n ) the operators used are 

A = Ai-i, Q — Qa, i» D = Da,u P — P A,i BJid R = Ra ,»• 

Since A n is diagonal, the method is well-defined. The total work of MBF for a problem with iV 2 
variables is 

w{N 2 ) = 0(N 2 ) + w(N 2 / 4) = 0(N 2 ) + 0(lV 2 /4) + tn(AT 2 /16) = • • • = 0(N 2 ) 

Note that if we had 

(rowsum(Rr,i+ 1 ) • rowsum(PT,i+ 1 )) 0 1 = Jo J) 5 ( j +1 
1 0 (rowsum(R s ,i+i) • rowsum(P s ,i+i)) = 0 I 


573 


then 


Qa,i+i = ^r,i+i2j+i 0 Es + Et 0 Ds,i+iSi+\ 

= Tj+i o Es,j + Ej,i o 5j+i 

is already in the scaled form, on which the MBF algorithm may act recursively. Hence the scaling 
by D~ l in action (4) of Section 2.2 is not needed. Of course, these equalities cannot hold exactly, but 
if they hold approximately, we can avoid scaling. Especially in the non-separable case, where scaling 
is impossible, action (4) has to take place without scaling by D~ l (instead of actual scaling, we would 
prefer in that case to keep diagonal matrices multiplying the difference operators in each of the space 
directions). Hence we would like to assume the above equalities at least for all former levels, that is, 
at the i ih call to MBF , for all 0 < j < i — 2. We call that variant K noscal . 

If we use the rectangular matrices P and R of the last part of Section 2.3, we get a variant of 
MBF which we call variant notri. In this variant, only the last of the four subsystems described 
above is solved. It assumes that the other tridiagonal subsystems affect only the smoothness of the 
error. The row-sum operation, however, is still performed on the original triangular matrices and not 
on the newly defined rectangular matrices. 

Instead of the operators and Q A ,i defined above, one may use difference operators which arise 
from the original PDE. If the algorithm is to be automatic, all such operators have to be of the 
same type (i.e. central or upwind) as the original fine- grid operator. (Nevertheless, in Section 3 
we will see that for some non-symmetric problems this condition has to be violated for the sake of 
convergence.) We use this strategy with the rectangular matrices of variant notri\ our version is then 
different from classical multigrid only in the choice of restriction and prolongation operators. Note 
that the row-sums computed in the MBF algorithm are usually 4. Instead of the multiplication by 
these row-sums, one may divide the residual by 4 before action (4) of Section 2.2. Then one gets an 
algorithm which is equivalent to that of [12] for the Poisson equation, and is close to that of [5] for 
other problems. We denote that strategy MGF (Multigrid + Factorization). It should be kept in 
mind that when applying this strategy one must use 2 n - 1 grid points on the finest grid and 2^ — 1, 
1 < q < fi for coarser grids in order to conserve uniformity. Here the even points, which are taken as 
coarse grid points, are always internal points of the original stencil. For 2 q point grids, on the other 
hand, the last fine grid point appears as a last grid point in all grids. Hence, coarse grids are biased 
towards the boundary. For our method MBF , on the other hand, stencils of both 2 n points or 2 n — 1 
points may be used. This is critical for implementation to problems on general regions, where grid 
lines may contain variable numbers of grid points (see Section 4). 

As mentioned above, the MGF method requires division of residuals by 4 before action (4) takes 
place. Sometimes It" is better to scale the discrete operators on all grids instead of dividing the 
residuals by this factor. Actually, for the Poisson equation both manners are equivalent: suppose A\ 
has a coarse step-size H = 2h ; then normalizing Ai to have the same diagonal entries as A$ = A 
amounts to multiplication of A\ by the factor 4, which is equivalent to the division of the residual 
in action (4) by that factor. Nevertheless, for differential equations that include derivatives of orders 
other than 2, this variant is not equivalent to MGF . We call it MGN (Multigrid + Normalization). 

The generalization of the MBF method and of the other multigrid versions to nonseparable 
problems is straightforward. A tensor product by an N x N diagonal matrix is to be replaced with 
a multiplication by an N 2 X N 2 diagonal matrix. 

Another generalization of MBF is to non-rect angular domains. This is also straightforward, since 
a line containing an odd number of points may be divided into two sets, one containing odd points 
and the other containing even points. Then one of those sets is considered as a coarse grid, and is 
divided recursively in the same way. A similar strategy may be used in the other space direction. 
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3 NUMERICAL EXPERIMENTS 


In this section, the MBF method is compared to other multigrid versions. The problems solved 
are of the type 

Lu{x,y) = f(x,y) (x,y) € (0,1) 2 
u(x,y) = xy (x,y) 6 d[0,l] 2 

The equations are discretized via a 3-point central difference scheme. For MBF, the number of grid 
points in each space direction is N = 2 n . For other multigrid strategies, however, the number of 
grid points in each space direction is N — 1 = 2” — 1; otherwise, coarse grids are biased towards the 
boundary (see Section 2.4), and the convergence is slow. 

For MBF, we have used variant “ notri ” of Section 2.4, in which tridiagonal subsystems are not 
solved. The main variant, which involves the solution of those subsystems, was found to be at most 
as effective as “notri”. 

The main variant of Section 2.4 was used for the hyperbolic, the non-uniform, the strongly in- 
definite and the discontinuous problems (the latter when applied with a red-black smoother). For 
other problems we have found that variant “noscal” described there (in which scaling of coarse-grid 
operators is omited), performs equally well. Hence we have chosen to use this simpler version rather 
than the main variant. Indeed, it was found that for most problems its performance was very close 
to that of the main variant. For the hyperbolic problem, however, its performance was twice as slow. 

The smoother of the error in all grids was the one provided by the 7X17(1,1) iteration of [24] 
[25] [26] or the red-black (RB) iteration. This determines the operators Z\ of Section 2.2 to be the 
preconaitioners for the 7X7(1, 1) or RB iteration, respectively. These smoothers were found to be 
superior to the Jacobi and damped- Jacobi smoothers. One presmoothing and one post-smoothing is 
performed. The initial guess is random. Double precision arithmetic is used. 

The integers in the following tables present the number of iterations needed to reduce the 1 2 norm 
of the residual by 10 6 . The maximum norm of the error was also computed, and its rate of convergence 
was close to that of the residual. 

In conjunction with the MBF iteration, we have used the computer code of [23] that implements 
the vector acceleration RRE that was mentioned in the introduction. The RRE acceleration was 
employed in cycling mode, by restarting it after every 10 iterations until convergence. The results of 
this are compared to those provided by the MBF iteration without acceleration denoted by NONE. 

We have also examined the classical multigrid versions mentioned at the end of Section 2.4. This 
strategy is denoted by the superscript D . The number of grid points in each space direction is 
TV — 1 = 2" — 1 . In most of the problems, the MGF and MGN versions of Section 2.4 are equivalent. 
Where this is not the case, we mention explicitly which of the two versions has been employed. 

For comparison, we also checked the performance of a method which does not involve any multigrid 
strategy. This method is the Modified 7X7 method of [29] with the optimal parameter of [30], used 
as a preconditioner for RRE. It is denoted by MILU. 

In the following tables, when a method converges very slowly we denote it by “slow”, and when 
a method diverges, we denote it by “div”. 


3.1 The Poisson Equation 


In table 1 we present the results for the Poisson equation. 
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RRE 

NONE 

RRE 

NONE 

RRE 

NONE 

RRE 

NONE 

RRE 

N 

ILU 

ILU 

RB 

RB 

mp 

ILU U 

RB U 

RB U 

MILU 

32 

4 

5 

5 

7 

4 

5 

5 

7 

19 

64 

4 

5 

5 

7 

4 

5 

5 

7 

27 

128 

4 

5 

5 

7 

4 

5 

5 

7 

42 


Table 1: Results for the Poisson equation 


MBF and the classical multigrid version perform equally well for this problem. 


3.2 Poisson Equation on a TchebychefF-Type Grid 


In table 2 we present the results for the Poisson equation, discretized via central differences on 
the 2-dimensional grid 


P 0 '.*) 3 ( 




1 < j,k < N 


The matrix operator for this scheme may be used as a preconditioner for a Tchebycheff-collocation 
discretization of the Poisson equation (see [31] and the references therein). 



RRE 

NONE 

RRE 

NONE 

RRE 

NONE 

RRE 

NONE 


N 

ILU 

ILU 

RB 

RB 


msam 

m wsm 

mmum 

MILU 

32 

4 

5 

7 

10 

4 

5 

8 

ii 

16 

64 

4 

6 

9 

19 

5 

6 

10 

18 

23 

128 

5 

7 

13 

33 

5 

6 

14 

28 

36 


Table 2: The Poisson equation with non-uniform grid 


The superscript D refers to the MGF method of the end of Section 2.4, which is in the spirit of 
Dendy [12J. It performs equally well as MBF. 


3.3 An Anisotropic Discontinuous Equation 


In table 3 we present the results for an anisotropic equation whose coefficients are discontinuous; 

a ( x ) u «* + o{y)u vy = 0 


Here a(t ) is defined by 



0.01 0 < t < 0.5 
1 0.5 < < < 1 


MBF and MGF perform equally well for this problem. For both methods, N — 1 grid points were 
used in each space direction. If iV = 2 n grid points are used in any of the space directions, then 
for coarse-grid problems the discontinuity lines are biased towards the boundary, and convergence 
becomes slow. 


Results similar to those of table 3 were obtained for the continues anisotropic problem 

Ujpjp “I” O.OlUyy — 0 

This time, however, there was no difference in convergence rate when the number of grid points in 
each space direction was changed from N = 2 n to N — 1 = 2" — 1 , for both MBF and MGF. 
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RRE 

NONE 

RRE 

NONE 

RRE 

NONE 

RRE 

NONE 

RRE 

N 

ILU 

ILU 

RB 

RB 

FWBS1 

mmm 

RB U 

RB U 

MILU 

32 

5 

6 

19 

slow 

5 

6 

19 

slow 

23 

64 

7 

10 

26 

slow 

7 

10 

26 

slow 

32 

128 

9 

13 

28 

slow 

9 

13 

28 

slow 

48 


Table 3: An anisotropic discontinuous equation 


3.4 A Convection-Diffusion Equation with Circular Streamlines 

In table 4 we present the results for the convection-diffusion equation 

u xx + Uy V + 150( (y — 0.5)u x — (x — 0.5)u v ) = / 
whose characteristics are circles: 



RRE 

NONE 

RRE 

NONE 

RRE 

NONE 

RRE 

NONE 

RRE 

N 

ILU 

ILU 

RB 

RB 


ILU U 

RB U 

RB° 

MILU 

32 

7 

10 

12 

15 

8 

slow 

div 

div 

28 

64 

6 

10 

9 

12 

8 

slow 

div 

div 

51 

128 

6 

10 

9 

12 

8 

slow 

div 

div 

94 


Table 4: A convection-diffusion equation with circular streamlines 


Problems of the last type are widely discussed in [6]. The approach developed there requires 
special treatments and is not as automatic as ours. 


3.5 A Convection-Diffusion Equation with Radial Streamlines 


In table 5 we present the results for the convection-diffusion equation 

u„ + Uyy + 150( xu x + yu y ) = / 
whose characteristics are radial lines: 



RRE 

NONE 

RRE 

NONE 

RRE 

NONE 

RRE 

NONE 

RRE 

N 

ILU 

ILU 

RB 

RB 

~IUJ V ~ 

mP 

RB U 

RB U 

MILU 

32 

7 

9 

slow 

div 

7 

9 

div 

div 

28 

64 

5 

6 

13 

12 

7 

8 

15 

15 

51 

128 

5 

6 

9 

10 

8 

9 

12 

15 

94 


Table 5: A convection-diffusion equation with radial streamlines 


The d superscript denotes here the A IGF method of the end of Section 2.4. Nevertheless, the 
purely automatic MGF version, in which all difference operators are central, diverged. To avoid 
that we had to use upwind difference schemes for all grids coarser than the original grid. This 
strategy, however, though performing almost equally well as MBF , suffers the disadvantage of not 
being automatic. Another way to overcome divergence is to use the MGN method of the end of 
Section 2.4, with the same step-size h for all grids. This strategy is non-automatic as well, and about 
twice as slow as the first one. 
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3.6 A Convection Equation 

In table 6 we present the results for the convection equation 

(y - 0.5 )u* - (x - 0.5)u y = f 

whose characteristics are circles, discretized via an upwind scheme: 



RRE 

NONE 

RRE 

NONE 

N 

ILU 

ILU 

ILU U 

TUP 

32 

10 

27 

10 

slow 

64 

15 

49 

14 

slow 

128 

23 

89 

23 

slow 


Table 6: A convection equation with circular streamlines 


The superscript D refers to the MGF method of the end of Section 2.4. Its performance is equal 
to that of MBF. 

The standard ILU and the Modified ILU of [29] (on the fine grid only, without multigrid strategy) 
do not converge for this problem. All multigrid strategies with an RB smoother are very slow. 


3.7 The Helmholtz Equation 


In table 7 we present the results for the Helmholtz equation 

U XX + Uyy + flu = f 

with (5 = 64. The RRE method for MBF was restarted in this example after every 5 iterations. 



RRE 

RRE 

RRE 

RRE 

N 

ILU 

RB 

TUP 

TOP 

32 

9 

14 

17 

16 

64 

9 

13 

17 

19 

128 

8 

15 

18 

20 


Table 7 : The Helmholtz equation 

Without acceleration, all methods diverged. The RRE acceleration for ILU and MILU iterations 
(on the fine grid only, without multigrid strategy) also diverged. 

The superscript D denotes here the MGN version, used with a continuation strategy; that is, use 
a parameter /3 smaller than that of the original PDE for grids coarser than the original grid, in such 
a way that the number h 2 Q is constant for all grids. Without this continuation strategy, divergence 
was reported. Consequently, it suffers the disadvantage of not being automatic. 

This problem is of the type of problems discussed in [7]. The projection approach given there 
requires more work and special treatment. 

For j3 > 64, the RRE acceleration for MBF seems to suffer stability problems, as the residual 
no longer decreases monotonically. A machine with higher precision (or Kacmarz smoother as in 
Section 3.8) is required. With classical multigrid, on the other hand, the acceleration is more stable: 
with the ILU smoother, it converges for /3 = 100 and N = 64 in 42 iterations. 

Note that the Helmholtz equation and the convection-diffusion equations are better-posed as the 
number of grid points increase; hence the number of iterations generally decrease. 
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3.8 Helmholtz Equation with Mixed Boundary Conditions 

The above experiments involve Dirichlet boundary conditions. In this sub-section we examine the 
Helmholtz equation with Dirichlet boundary conditions on the edges x - 0, x - 1 and y - i, and 
with the mixed boundary conditions 

+ = S (5) 

on 


on the edge y = 0. 

We have repeated experiments 6.1 and 6.3 of [32], for which /? = 100, a = lOOi, N = 31 and 
8 = 200 a = 10 i, N = 63 respectively. On coarse grids, where the problem is very indefinite, we 
have used Kacmarz relaxations as a smoother. The cost of an MBF or MGF iteration was about 
10 Jacobi iterations of the original problem. With MBF , we have converged in 10 iterations for the 
first problem and in 18 for the second one, which is much better than the results of [32]. With MG* , 
applied with the continuation strategy described in Section 3.7, we have converged m 23 iterations 

for the first problem. 


3.9 Helmholtz Equation with Finite Elements 


Finally, we have examined the Helmholtz equation 


u 3 


+ Uyy + fa = / 


with 8 = 200 and Dirichlet boundary conditions, discretized via bilinear finite elements. The grid 
for those elements is not uniform; in each space direction, the domain is divided into 4 elements, and 
the grid points are the roots of the Legendre polynomial of degree 17 in each element. Hence the 
total number of grid points is 63 2 . This grid induces a division of the domain into squares, which 
serve as the bilinear finite elements. The matrix operator for this bilinear element scheme may be 
used as a preconditioner for a spectral element discretization of the Helmholtz equation (see [31] and 
the references therein). Though the coefficient matrix has nine non-zero diagonals, the operators 
for coarser grids have five non-zero diagonals only; they are obtained from the above finite difference 
approximation in the automatic or classical manners. Actually, this is a double discretization strategy. 
The relaxations on the finest grid are the ILU iteration or the four-color Gauss-Seidel iteration. On 
the second grid, the relaxation is ILU or RB iteration. One presmoothing and one post-smoothing 
are performed on those two grids. On coarser grids, since the operators are more indefinite, these 
relaxation methods are too divergent; hence, we use instead the Kacmarz iteration 40 presmoothmgs 
and 40 post-smoothings on each level. Since on the third grid the number of points is 1/16 of that of 
the original grid, the total work on that grid is about five Jacobi relaxations of the original system. 
The cost of the whole multigrid or MBF procedure is about 10 such relaxations. RRE acceleration 
is restarted after every 10 multigrid or MBF iterations. The number of MBF iterations needed to 
reduce the residual by 6 orders of magnitude is 28 when ILU is used on the two finest grids and 27 
when the Gauss Seidel smoother is used there. For MGF (with ILU on the two finest grids and wi h 
the continuation strategy of Section 3.7), the number of iterations needed is 52. When the residual is 
reduced by 6 orders of magnitude, the error is reduced by 6 orders for MBF and 5 orders for MG* . 

We also examined the mixed boundary conditions case. For the mixed boundary conditions (5) on 
the edges x = 0 and y = 0, with 0 = 200, a = 10* and N = 64, MBF converged in 52 iterations, each 
costs about as much as 7 Jacobi iterations, with RRE restarted after every 20 iterations. Classical 
MG methods did not converge for this problem, even with the continuation strategy of bection o. (. 
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3.10 Helmholtz Equation with Spectral Elements 

r,J^±° V % d T h \ e d i s , cretizat i ion 1 st , rate «y is , not limited to bilinear element schemes; it was em- 
{ ? T ^Vu PeCtra elem ® nt scheme of [33] as well. Again, we relax the original 
quation on the finest grid, then we compute the difference scheme based on the nodes of the spectral 
dements, and use it to generate the coarse grid operators via MBF. These operators are used to 
J™ d ( he Jf arse g nd correction. For the Helmholtz equation with the mixed boundary conditions (5) 
on the edges x = 0 and y = 0, with the parameters 0 = 100, a = 100z' and N = 16 and with 

s P e ? t r a l e lem e nts, MBF converged in 16 iterations, each costs about as much 
as 5 Jacobi iterations (with RRE restarted after 10 iterations). For the same problem but with the 

I- 2 -? 0, \ = 10 Ti and Z 64 > MBF converged in 60 iterations; each costs about as 
much as 2 Jacobi iterations. This rate of convergence was much better than that of the algorithm of 
IdlJ, in which the spectral element scheme is preconditioned by a finite element or a finite difference 
scheme (even when the preconditioning system is solved by MBF). As a matter of fact, even for the 
Poisson equation our strategy was about 3 times faster than that of [31] 


4 DISCUSSION 

on met5 i od 1S a ve rsion of multigrid, which is automatic in the sense that it depends 

on the algebraic system of equations rather than on the original PDE. Actually, it is a “black box” 
method for the solution of the linear system of equations. Hence it seems to be more robust than the 
classical multigrid method. For instance, nonsymmetric terms in the equation do not slow down the 
onvergen.ee, whether the characteristics are closed or open. Non-uniform grids are handled with the 
same efficiency, and no special treatment of the neighborhood of the boundary is needed. Moreover 
^^ acceleration is applied to the method, it copes with the indefinite Helmholtz equation 
as well. For all these examples the rate of convergence is rather independent of the size of the problem 
ror anisotropic or P u I e advection problems, however, the rate of convergence of the MBF method 
applied with the RRE acceleration slightly depends on the size of the problem. 

The MBF method is especially suited to use with the ILU smoother. The red-black smoother 
gives slightly worse results. The doubled damped Jacobi iteration as a smoother (with a damping 
parameter 0.5) was examined too. For all the above problems but the discontinuous-anisotropicand 
the hyperbolic problems, its performance was about twice slower than that of the ILU smoother. For 
those two problems, the damped Jacobi smoother was unsatisfactory. 

The versions of multigrid denoted MGF and MGN perform well for problems which do not involve 
central first derivatives (including discontinuous and anisotropic problems). For problems which do 
“f, ta] £ ^ n .V al first deriv atives, since the algorithm is assumed automatic, the discretization on coarse 
^. e same type ^ that ° f the fine grid, i.e. central. Hence divergence is often caused by 

1 rs ®'g nds corrections. This difficulty can be handled by the special treatments of Section 3.5 

! he al .g° nthm » no onger automatic For the Helmholtz equation, one may overcome this 
difficulty by using a continuation strategy in the MGN version. Even though this (non-automatic) 
strategy is a bit slower than MBF , it is more stable and is applicable to more singular problems If 

b0th MBF and MGF "Sr to"** 

As opposed to the classical multigrid versions, the MBF is applicable whether the number of grid 
Ch - PaCe i rectl ° n 18 < r ven ° r odd o This indicates that it is applicable to problems defined 

2 dbn^tinaTfi 0118 ' p glon C R » .? ne take ? as a fine g rid the restriction of an infinite 

j' dimensional fine gnd to H. For a coarser grid, one takes every other point (in both x and y space 

directions) in the infinite fine grid, and takes the restriction to Cl. The other coarse grids are created 

nnw!t! ? ame W£ 2o A i We haVC Se ®j’ MBF 18 not affected b y the possibility that some coarse grid 
' VT ^ | coarse grid operators are created automatically as in the above description; 

tins can be done easily by modifying the block sizes in the coefficient matrix of the system. Thus the 
algorithm is easy to program. 7 3 lue 
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